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Abstract 

Kinetic properties of a granular gas of viscoelastic particles in a homoge- 
neous cooling state are studied analytically and numerically. We employ the 
most recent expression for the velocity-dependent restitution coefficient for 
colliding viscoelastic particles, which allows to describe systems with large 
inelasticity. In contrast to previous studies, the third coefficient a 3 of the 
Sonine polynomials expansion of the velocity distribution function is taken 
into account. We observe a complicated evolution of this coefficient. More- 
over, we find that 03 is always of the same order of magnitude as the leading 
second Sonine coefficient a 2 ; this contradicts the existing hypothesis that 
the subsequent Sonine coefficients 02, 03 . . ., are of an ascending order of a 
small parameter, characterizing particles inelasticity. We analyze evolution 
of the high-energy tail of the velocity distribution function. In particular, 
we study the time dependence of the tail amplitude and of the threshold 
velocity, which demarcates the main part of the velocity distribution and the 
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high-energy part. We also study evolution of the self-diffusion coefficient D 
and explore the impact of the third Sonine coefficient on the self-diffusion. 
Our analytical predictions for the third Sonine coefficient, threshold velocity 
and the self-diffusion coefficient are in a good agreement with the numerical 
finding. 

Key words: Granular gas; Viscoelastic particles; Velocity dependent 
coefficient of restitution; Velocity distribution function; Self-diffusion 



1. Introduction 

An ensemble of macroscopic particles, which move ballistically between 
dissipative collisions, is usually termed as a granular gas, e.g. 
and the loss of energy at impacts is quantified by a restitution coefficient e, 

(#12 ■ e) 



(1) 



(vi2 ■ e) 

Here V\2 and v[ 2 are relative velocities of two particles before and after an 
impact and e is the unit vector, connecting their centers at the collision 
instant. In what follows we assume that the particles are smooth spheres of 
diameter a, that is, we do not consider the tangential motion of particles. 
Then the restitution coefficient yields the after-collision velocities of particles 
in terms of the pre-collision ones, 



iTi - -(I +e) {V12 ■ e) e. 
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[1 + e) (v 12 -e)e. 



(2) 
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The simplest model for the restitution coefficient, e = const, facilitates signif- 
icantly the theoretical analysis and often leads to qualitatively valid results, 
e.g. y). This assumption, however, contradicts the experimental observa- 
tions, e.g. 

Baa 

and basic mechanical laws, e.g. [9|, [10|. A first-principle 



analysis of a dissipative collision may be performed, leading to a conclu- 
sion that e must depend on the relative velocity of the colliding particles 
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12] . If the impact velocity is not very high (to avoid plastic deforma- 
tion of particles material) and not too small (to avoid adhesive interactions 



at a collision 13|), the viscoelastic contact model may be applied 11] . This 
yields the velocity-dependent restitution coefficient fl4. \\\: 

e = 1 - K6 (2n(t)) 1/10 |(c 12 e)| 1/5 + C 2 5 2 (2n(t)) 1/5 |(c 12 e)| 2/5 ± . . . (3) 

Here C x ~ 1.15, C 2 ^ 0.798 and u(t) = T(t)/T(0) is the dimensionless 
temperature, expressed in terms of current granular temperature, 

-nT{t) = J dv—f(v, t) = , (4) 

with f(v,t) being the velocity distribution function of grains, n is a gas 
number density, m is the particle mass and vr{t) is the thermal velocity; 
C12 = vi 2 /vt is, correspondingly, the dimensionless impact velocity. 

The restitution coefficient e ^ depends on the small dissipation param- 
eter 5 
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s^my, (5) 

\ m J 

proportional to the dissipative constant A, 

l (3 V 2-Vi) 2 (l-v 2 )(l-2v) 
3 (37/2 + 277O Yv 2 

which depends on the Young modulus of the particle's material Y, its Poisson 

n 

ratio v and the viscous constants r\\ and 772, see e.g. [1|. The parameter k 



Due to inelastic nature of inter-particle collisions, the velocity distribution 
function of grains deviates from the Maxwellian, so that the dimensionless 
distribution function ill 



71 ~ 

f(v,t) = ^f(c,t) 



may be represented in terms of the Sonine polynomials expansion, e.g. 



(8) 
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221: 




(9) 



/(c,t) = 0(c) \l + ^2a p (t)S. 

P =i 

where d>(c) = 7r~ 3//2 exp(— c 2 ) is the dimensionless Maxwellian distribution 

Q. 



and the first few Sonine polynomials read 



S (x) 
S 2 (x) 



1, S 1 (x) 
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Provided that the expansion (J5]) converges, the Sonine coefficients a p com- 
pletely determine the form of f(c,t). According to the definition of temper- 
ature ai(t) — 1|, that is, the first non-trivial coefficient in the expansion 
is a,2{t). For the case of a constant e the coefficient a>2 has been found 




v | X 61 . 117IJ and the coefficient 03 analytically [18( and numerically 
2l| . The expansion (Q quantifies the deviation of the distribution 
function from the Maxwellian for the main part of the velocity distribution, 
that is, for c ~ 1; the hi gh-v elocity tail of f(c) for c ^> 1 is exponentially 
over-populated , I22J and requires a separate analysis. 



The impact- velocity dependence of the restitution coefficient, as it fol- 
lows from the realistic visco-elastic collision model, has a drastic impact on 



the granular gas properties. Namely, the form of the 



function and its time dependence significantly change 
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velocity distribution 



23], similarly changes 



251 ] . Moreover, the global 



the time dependence of the kinetic coefficients 
behavior of the system qualitatively alters: Instead of evolving to a highly 
non-uniform final state of a rarified gas and dense clusters, as predicted for 



e = const. [26|, |27|, the clustering [27j and the vortex formation 



28 



occurs 



in a gas of viscoelastic particles only as a transient phenomena [29|. 

Based on the restitution coefficient given by Eq. ([3]) the theory of granular 
gases of viscoelastic particles has been developed, e.g. [lj], where as in the 
case of a constant e, only the first non-trivial Sonine coefficient a-i (t) has been 
taken into account. Although the evolution of the high-velocity tail of f(c) 



231 ] . neither the location of the tail, nor its amplitude was 



has been analyzed 
quantified. 

Recently, however, it has been shown for the case of a constant restitu- 
tion coefficient, that the next Sonine coefficient 03 is of the same order of 
magnitude as the main coefficient 181 ]: moreover, it was also shown that 



the amplitude of the high-velocity tail of f(c) and its contribution to the 
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3l|- 



kinetic coefficients may be described quantitatively 

Finally, a new expression for the velocity-dependent restitution coefficient 
e has been derived 32j . which takes into account the effect of "delayed re- 



covery" in a collision. The delayed recovery implies, that at the very end 
of an impact, when the colliding particles have already lost their contact, 



their material remains deformed 



321 ] . This affects the total dissipation at a 
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collision, so that the revised e reads 



32]: 



e = l + J2^ /2 (2u(t)y /20 \(c l2 e)\ l/w . (13) 
i=i 

Here hi = 0, h 2 = —C\, h 3 — 0, — C 2 and the other numerical coefficients 
up to i = 20 are given in |32j. Fig. [1] (inset) illustrates the dependence of 



the restitution coefficient (fT3|) on the dissipative parameter S for u = 1 for 
collisions with the characteristic thermal velocity, that is for |(c 12 • e) \ = 1. 

In the present study we use the revised expression ([TBI for e and investi- 
gate the evolution of the velocity distribution function in a gas of viscoelastic 
particles in a homogeneous cooling state. With the new restitution coefficient 
we are able to describe collisions with significantly larger dissipation, than 
it was possible before with the previously available expression for e. For the 
larger inelasticity one expects the increasing importance of the next-order 
terms in the Sonine polynomials expansion and of the high-energy tail of the 
velocity distribution. In what follows we study analytically and numerically 
time evolution of f(c,t), using the Sonine expansion up to the third-order 
term and analyze the amplitude and slope of the high-energy tail of f(c,t). 
In addition, we consider self-diffusion - the only non-trivial transport process 
in the homogeneous cooling state and compute the respective coefficient. 

The rest of the paper is organized as follows. In the next Sec. II we 
address the Sonine polynomial expansion and calculate the time-dependent 
Sonine coefficients, along with the granular temperature. In Sec. Ill the 
high-energy tail is analyzed, while Sec. IV is devoted to the self-diffusion 
coefficient. Finally, in Sec. V we summarize our findings. 
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2. Sonine polynomial expansion: evolution of the expansion coef- 
ficients 

Evolution of the velocity distribution function of a granular gas of spheri- 
cal particles of diameter a in a homogeneous cooling state obeys the Boltzmann- 
Enskog equation, e.g. 

^|^ = »M/(/,/), (14) 

where /(/, /) is the collision integral. Generally, I depends on the two- 
particle distribution function ^(^i, v 2 , r 12 , t). Within the hypothesis of molec- 
ular chaos, f 2 (vi,v 2 ,f 12 ,t) = g 2 (a)f(vx,t)f(v 2 ,t) (see e.g. [lj), the closed- 
form equation (I14p for f(v, t) is obtained; g 2 (a) denotes here a pair correlation 
function at a contact, which accounts for the increasing collision frequency 
due to the effect of excluded volume [1]. 

Eq. (THI) yields the equation for the dimensionless distribution function 



f(c, t) [l| 



y (3 + cj£\ f(S, t) + B'^f(c, t) = /(/, /) , (15) 



where B = vrg2(cr)a 2 n and 

th = ~ J d<&i{f,f), (16) 
is the p-th moment of the dimensionless collision integral: 

/(/,/) = j dc 2 J deQ(-c 12 -e) |-c 12 ■ e| ^/(c/', t)/(c 2 ", t) - f(c u t)f(c 2 ,t)^ 

(17) 

Q(x) in the above equation is the Heaviside step-function and the dimen- 
sionless velocities c" and c 2 are the pre-collision velocities in the so-called 



inverse collision, which results with c\ and c 2 as the after-collision velocities, 
e.g. fl. 

Eq. Q15p is coupled to the equation for the granular temperature, e.g. [l| : 



dT 2 



(18) 



which defines the cooling coefficient ( = (2/3)Bfi 2 - 

Multiplying Eq. ffT5]) with c p , integrating over c and writing Eq. (TT8]) in 
the dimensionless form, we obtain for p — 4, 6 jlj] (similar equations has been 



used in 



du 
dr 

da? 
dr 

das, 
dr 



^/i 2 (l + a 2 )- T ^/i 4v ^ 
^ 2 (1 - a 2 + a 8 ) - ^/i 4 v^ + 



(19) 



Here r = t/r c (0) is the dimensionless time, measured in the mean collision 



units, r~ 1 (0) = A^g2(a)a 2 n^T(0)/m at initial time t = 0. The equations 
for the Sonine coefficients a 2 and 03 in ([19 ) are the first two of the infinite 
system of equations for p = 4, 6, 8, . . . [34j, QJ. The moments /z 2 , /i4, /i6, ... in 
these equations depend on all the Sonine coefficients a^, with k = 2, . . . , 00, 
that is, only the infinite set of equation is closed. To make the system 
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33[ we truncate the 



tractable one has to truncate it. Following Refs. 
Sonine series at the third term, i.e. we approximate, 

/(c, r) ~ 0(c) (1 + a 2 (r)S 2 (c 2 ) + a 3 (r)S 3 (c 2 )) . (20) 

Within this approximation /i 2 , /i 4 , /x 6 in Eqs. ( 1T91) depend only on a 2 , 03 and 
n, which makes the system ffl9l) closed. 

The moments of the collisional integral fi p were calculated analytically up 
to 0(5 10 ), using the formula manipulation program, explained in detail in 
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The complete expressions are rather cumbersome, therefore, we present 
here for illustration only linear approximations with respect to a 2 , a 3 and 
6'(t) = 5(2u(r)) 1/10 : 



6 2 
H2 = u [l + —a 2 + — a 3 ) 5\t) 



/4 129 179 \ 
fi4, = V27T (4a 2 - a 3 ) + 7 ( - + — a 2 - — a 3 J cV(r) (21) 

/^-/ ir 45 \ , /1079 40717 39353 \ 

«6 = 3v 2?t 15a 2 a 3 + 3cj o (t) 1 a 2 a 3 , 

Pb V 4 7 V 100 1250 3125 J 

where w = 2v / 27r2 1 / 10 r (21/10) d ~ 6.485. 

Since the system of equations ( fl9l) is strongly non-linear, only a numerical 

solution is possible. Still, one can find a perturbative solution in terms of 

small dissipative parameter 5, 

a 2 = a 20 + a 2 ±S + a 3 = a 30 + a 3 iS + u = u + uiS+... (22) 

The zeroth-order solution (with a 2 (0) = a 3 (0) = 0) reads, 

a 20 = 0, ^30 = 0, m = (1 + r/T y 5/3 , (23) 

where Tq 1 = 2 6 / 5 r (21/10) C\5/h ~ 0.555. For the first-order solution only 
its asymptotics for r — ► 00 may be obtained analytically: 

1 57 

a 2 = -A 2 5(r/r )- 1/6 A 2 = 2 1/5 ^T (21/10) C x ~ 0.44 (24) 

a 3 = -A 3 5(r/r )^ 1/6 A 3 = 2 1/5 ^r (21/10) C x ~ 0.08 (25) 
^ = (r/ro)^ 5/3 + g5(r/r )- 11/6 , (26) 
^2^^1(21/10) + 1^^3.28. 

y V 15625 r (21/10)7 

Interestingly, the third Sonine coefficient a 3 is of the same order of magnitude 
with respect to the small parameter 5, as the second Sonine coefficient a 2 , 



albeit five times smaller. This conclusion is in a sharp contrast with the 
hypothesis suggested in Ref. [34|, that the Sonine coefficients are of an 
ascending order of some small parameter A, that is, ~ A . 




Figure 1: Evolution of the dimensionless granular temperature u(t) = T(r)/T(0), as 
it follows from the numerical solution of Eqs. (fT9|) . for different values of the dissipativc 
parameter 5. Time is measured in the collision units tq (see the text). For t 3> 1 all curves 
demonstrate the same slope, u(t) ~ t -5 / 3 (shown by the dotted line), in accordance with 
the theoretical prediction, Eq. ([26]) . The inset illustrates the dependence (fT3|) of the 
restitution coefficient e on the dissipativc parameter 5 for u = 1 for collisions with the 
characteristic thermal velocity, that is for \c\2 ■ e \ = 1. The particular value of e = 0.7, 
discussed below, corresponds to <5 ~ 0.3. 

The numerical solution of Eqs. (jT9{) confirms the obtained asymptotic de- 
pendence, Eqs. (!24l) - ([26l) . This is seen in Fig. [fl where the time dependence 
of the reduced temperature u{r) is plotted and in Figs. [2] and [3] (insets). Fig. 
CD demonstrates that the larger the dissipative parameter 5, the earlier the 
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asymptotic behavior of u(t) is achieved. 

The numerical solution for the Sonine coefficients 02 (t) and 03 (r), shown 
in Fig. [2] and [31 corresponds to the initial Maxwellian distribution, 02(0) = 
°3(0) = 0. As it follows from the figures, the absolute values of the both 
coefficients initially increase, reach their maxima |a2max| and |ct3 max | and, 
eventually, decay to zero. In other words, the velocity distribution function 
for viscoelastic particles evolves towards the Maxwellian. It is interesting to 
note, that the maxima |a2 max | and 1 0,3 mELX | first increase with the increasing 
dissipative parameter 5, then saturate at 5 ~ 0.15 and do not anymore 
grow. The location of the maxima, however, shifts with increasing 5 to later 
time, Figs. [2] and [31 This illustrates the general tendency - the larger the 
dissipation parameter, the slower the gas evolution. Again we see that the 
third Sonine coefficient (13 is of the same order of magnitude as 02 (although 
a few times smaller) for all evolution stages. 

To understand the observed behavior of the Sonine coefficients, consider 
the dependence of these coefficients on e for a constant restitution coefficient 
(Fig. HJ) and the respective dependence of the restitution coefficient e on the 
dissipative parameter 5 (Fig. [TJ inset). For < 5 < 0.15 [which corresponds 
to 0.85 < e < 1, Fig. [T] (inset)] one observes a fast relaxation, on a collision 
time scale, of the Sonine coefficients 02, 03 to their maximal values |a2max|, 
1 03 max I, roughly corresponding to the respective values for the constant e. In 
a course of time the granular temperature u(t) decreases and the effective 
restitution coefficient alters in accordance with decreasing dissipative param- 
eter S'(t) = 5 (2n(r)) 1 ^ 10 , that is, the effective e increases with time. For 
this interval of e (0.85 < e < 1) the increasing restitution coefficient implies 
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the decrease of the absolute values of the Sonine coefficients, Fig. HJ which 
is indeed observed in the evolution of a 2 (r) and 03 (r). The larger the 5 (i.e. 
the smaller the effective restitution coefficient), the larger the maxima |a2max| 
and |a 3max |, in accordance with the dependencies a 2 (e) and (13(e) depicted in 
Fig. U 

Similarly, for 0.15 < 5 < 0.30 [which corresponds to 0.7 < e < 0.85, Fig. 
[Tj (inset)] the initial fast relaxation to the related values of |a2max| and |a3 max | 
first takes place. Then the decreases of S'(t) and the respective increase of 
the effective e causes the increase of the absolute value of a 2 (t), until it 
reaches the maximum, corresponding to S'(r) = 0.15 (or e = 0.85), Fig. 
HI The further decrease of S'(t) leads to the according decay of a 2 (r), in 
agreement with the dependence of a 2 (r) shown in Fig. [2J This qualitatively 
explains the evolution of a 2 (r) and the saturation of its maximum |a2max| 
for S > 0.15. The qualitative behavior of the third Sonine coefficient may be 
explained analogously. 

It is interesting to note that the observed dependence of a 2 (r) with the 
new restitution coefficient ffTB"]) differs qualitatively for 5 > 0.145 from that 
obtained previously for the old restitution coefficient ([3]). While in the latter 
positive bump at initial time, r ~ 1 was detected for 5 > 0.145 



231 ]. in the former case the positive bump appears at much larger 5 > 0.3, 



Fig. [51 This is again in agreement with the behavior, expected from the 
dependence of a 2 (e) for a constant e: The coefficient a 2 becomes positive 
for e < 0.7, which corresponds to 5 > 0.3, Fig. [1] (inset). Note, however, 
that with increasing 5 more and more terms in the expansion ([TBI are to 
be kept. While for 5 = 0.27 it suffice to keep terms up to 5 9 , for 5 = 0.33 
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Figure 2: Evolution of second Sonine coefficient 02 (r). At the first stage of the evolution 
1 02(7") | increases, then reaches the maximum value |a2max|, and eventually decays to zero. 
With the increasing dissipation parameter S the maximum |a2max| shifts to later time. The 
inset shows the asymptotic dependence of 1 02(7") I ( m ll line) together with the analytical 
result (|24|) (dotted line) for S — 0.01 in the log-log scale. 
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Figure 3: Evolution of the third Sonine coefficient 0,3 (r). The inset shows the asymptotic 
dependence of |a3(r)| (full line) together with the analytical result (|25|) (dotted line) for 
S = 0.01 in the log-log scale. 
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Figure 4: Dependence of the second and third (inset) Sonine coefficients on the restitution 
coefficient e for the case of a constant e Q, 0, H]- Note that £ ~ 0.7 corresponds to the 
vanishing ct2- 

the Sonine coefficients a-i computed with the accuracy 0(5 9 ), 0(5 19 ^ 2 ) and 
0(5 10 ) noticeably differ, Fig. Therefore we conclude that the revised 
restitution coefficient for a viscoelastic impact (TT3T) may be accurately used 
up to 5 = 0.3. The loss of the accuracy in the computation of the Sonine 
coefficients for larger values of 6 may be a manifestation of the breakdown 
of the Sonine polynomials expansion, as it has been found previously for a 



constant restitution coefficient [18]. Similarly as in the case of constant e 
jisj], we expect that in the domain of convergence of the Sonine expansion, 
5 < 0.3, the magnitude of the next-order Sonine coefficients CI4, CI5, ... is very 
small, so that an acceptable accuracy may be achieved with the use of the 
two coefficients, 02 and 03 only. 
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Figure 5: The time dependence of the second Sonine coefficient a 2 (r), computed with a 
different number of terms in the expansion (fT5|) for the restitution coefficient. The full 
lines correspond to the accuracy up to 0(S 10 ), the dashed lines - up to 0(5 19 ^ 2 ) and the 
dotted lines - up to 0(S 9 ). Note that for S = 0.33 the accuracy up to 0(<5 10 ) is insufficient 
to obtain a reliable convergence. 
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3. High-velocity tail of the velocity distribution function 

The expansion ([9]) refers to the main part of the velocity distribution, 
c ~ 1, that is, to the velocities, close to the thermal one, vj. The high- 



velocity tail c ^> 1 is however exponentially overpopulatec 



in a course of time, during a first few tens of collisions 
particles the velocity distribution reads for c ^> 1 

/(c,r) ~ exp (-p(r)c) , 
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15| . It develops 



31j . For viscoelastic 



where the function ip(r) satisfies the equation 23| 



(27) 



:28) 



with B and /i 2 defined previously. Using ^ 2 (t), obtained by the formula 
manipulation program |l| (see the discussion above), we solve numerically 
Eq. (j28p to obtain ip(r). In the linear approximation, /i 2 « 6.49 5 (2m(t)) 1 ^ 10 , 
the function <p(r) has the form: 



<p(T) = (b/6)(l + T/T )* 



(29) 



with b w 1.129 ll|. 

Following [30j we neglect the transition region between the main part of 
the velocity distribution function, c ~ 1, and its high-energy part, c> 1 and 
write the distribution function as 

f(c,r) = A(r)c 2 exp(-c 2 ) (1 + a 2 (r)S 2 (c 2 ) + a 3 (r)5 3 (c 2 )) 6(c* - c) + 

+ J B(r)c 2 exp(-v?(r)c)e(c-c*) (30) 

The coefficients A(t), B(t) and the threshold velocity c*, which separates 
the main and the tail part of f(c, r) can be obtained, using the normalization 
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condition: 

J f(c)dc = 1 (31) 
and the continuity condition for the function itself and its first derivative 



214 : 



/( c *-0,r) = /(c* + J r) (32) 

0/(c*-O,r) _ 3/> + 0,t) 

& " dc {66) 

Substituting Eq. (JUD into (J5I|), (j32j) and (|53)l we arrive at, 

' (2c* -<p)(l + a 2 S 2 (c* 2 ) + a 3 S 3 (c* 2 )) = a 2 (2c* 3 - 5c*) + a 3 (7c* 3 - c* 5 - f c*) 
£ = Aexp(-c* 2 + ^c*) (1 + a 2 S 2 (c* 2 ) + a 3 S 3 (c* 2 )) 

(34) 

(12v^Ferf(c*) exp(c* 2 ) - 12a 2 c* 5 + 30a 2 c 3 * + 35a 3 c* 3 - 28a 3 c* 5 + Aa 3 c* 7 - 24c*) 
k x A eX p(-c* 2 ) + ^ (2 + <^c*(2 + <pc*)) exp(-y? C *) = 1 

where <p(t) is the solution of Eq. (|28p . 

To find the amplitudes A(r) and B(t) together with the threshold velocity 
c* the system (134)) was solved numerically. 

The asymptotic dependence of c* on r may be easily found if we take into 
account that a 2 and a 3 are of the same order of magnitude for r 3> 1, while 
<p(t) ^> 1 and c* 3> 1. Keeping in the first equation in (1341) only the largest 
terms, yields (2c* — </?)a 3 c* 6 ~ — a 3 c* 5 , which implies that 



c T 



^(r)/2 = (6/25)(l + r/r ) 1 / 6 (35) 



The typical velocity distribution function, computed for r = 50 and 5 = 
0.3 is shown in Fig. [61 The threshold velocity, c* separating the main and 
the tail part of the velocity distribution reads in this case, c* ~ 4.31. The 
threshold velocity increases with decreasing inelasticity 5 and shifts to larger 
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c 



Figure 6: A typical velocity distribution function in a granular gas of viscoclastic particles 
at r = 50 for S — 0.30 (Eq. (|5D|) . solid line), which is represented as a sum of velocity 
distribution function, obtained in the Sonine polynomial approximation (Eq. (|20p . dashed 
line), and the exponential function (Eq. (|27|) . dash-dotted line). For comparison the 
Maxwcllian distribution is also shown (dotted line). The threshold velocity, c* = 4.31, 
may be compared with the respective quantity c* = 3.77 for a gas of particles with a 
constant restitution coefficient e ~ 0.71. 
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Figure 7: Evolution of the threshold velocity c* in a gas of viscoelastic particles for S = 
0.3, 0.2, 0.15 (from top to bottom). The corresponding quantities for the case of a constant 
restitution coefficient arc: c*(e ~ 0.71) = 3.77, c*(e ~ 0.79) = 3.65 andc*(e ~ 0.84) = 3.60 
[30I ] - In a course of time the threshold velocity shifts to larger values, that is, the high- 
energy tail becomes less pronounced. The inset shows the asymptotic dependence of c* 
(full line) together with the analytical result ([33)) (dotted line) for 5 = 0.2 in the log-log 
scale. 
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values at later time, Fig. [3 which means that the high-energy tail becomes 
less pronounced. Again, we see that in contrast to a gas of particles with 
a constant restitution coefficient, where the tail of f(c) persists after its 
relaxation, in a gas of viscoelastic particles the velocity distribution function 
tends to a Maxwellian. 



4. Self-diffusion 

Self-diffusion is the only transport process which takes place in a granular 
gas in a homogeneous cooling state: In spite of the lack of macroscopic cur- 
rents, a current of tagged particles, identical to the particles of the surround- 
ing gas, but somehow marked, may exist. Moreover, the diffusion coefficient 
D is directly related to the mobility coefficient of the tagged particles k via 
the Einstein relation, k ~ D/T, which approximately holds true for granular 



361 ] . The mean-square displacement of the tagged particles 



gases, e.g. [35 
reads [l| 

<[Ar(r)] 2 > = ^D^df, (36) 

where the time-dependent diffusion coefficient (diffusivity) is the solution of 
the following equation [l|: 

-erf + ^-4. (37) 

Here £ = (2/3)Bfi 2 is, as previously, the cooling coefficient and T Vj0l d is 
the velocity correlation time, which characterizes the time after which the 
memory about the initial particle velocity is lost; it reads in terms of the 



distribution function 



1]: 



r v ,ad = ^v T {t)g 2 (a)a 2 n x (3* 
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Figure 8: Self-diffusion coefficient as a function of time, measured in the collision units 
To (see the text). The dissipative parameters, from the top to bottom are <5 = 0.3 and 
5 = 0.1. Dq is the Enskog self-diffusion coefficient in a gas of elastic particles at the initial 
temperature Tq = T(0). The full lines correspond to the complete solution of Eq. ([BT]) . 
while the dashed lines depict the linear approximation to D, (|42[) . Note that after a short 
period of time of about ten collisions, the complete and linear solution become practically 
indistinguishable. With the dotted lines we show u x ' 2 {t) = (T(t)/T ) 1/2 , which is equal 
to the ratio of two diffusion coefficients for elastic particles at temperatures T(r) and To- 
The plotted ratio D(t)/Do clearly tends to u 1 / 2 (r), which manifests that the self-diffusion 
in a gas of viscoelastic particles tends to that in a gas of elastic particles. The inset shows 
the relative deviation AD = (D — D)/D of the diffusion coefficient D computed with the 
use of both Sonine coefficients, ai and 03 from the respective value, D, obtained with the 
use of <i2 only (6 = 0.3). Naturally, the location of maximum of AD coincides with that 
of |o 3 (t)|, Fig.H 
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x J dcic 2 J de®(-c 12 ■ e) \c 12 ■ e\ f(c 1 ,t)f(c 2 ,t)(l - e)(c 12 ■ e) 2 . 

With the approximation ( 1201) for f(c,t) r v ^d was calculated up to 5 W , using 
the formula manipulation program, described in Here we present for 
illustration only its linear with respect to a 2 , and 5' part: 

( \ + A fl2 + L a \ _ (1+ ^-a 2 + ^-a 3 ) S'(t) (39) 



t e (0) \ 16 64 J \8 100 500 

Here rg(0) is the Enskog velocity correlation time of elastic particles at initial 
time r = 0: 

Using the obtained expressions for /i 2 and r Uia< i up to 0(5 10 ), we solve Eq. 
(I37p numerically and compute the diffusion coefficient D(t). In Fig. [H]the 
ratio of D(t)/D is plotted, where 

A> = r £ (0)T(0)/m, (41) 

is the Enskog self-diffusion coefficient for a gas of elastic particles at initial 
temperature T(0). As it is clearly seen from the figure, the diffusion co- 
efficient decreases with time in a way, that the ratio D{t)/Dq approaches 
= (T(r) /T(O)) 1 ^ 2 - the ratio of the two diffusion coefficients for elas- 
tic particles at the current temperature T(r) and the initial temperature 
T(0). Hence in a course of time the self-diffusion in a gas of viscoelastic 
particles tends to that in a gas of elastic particles. 

In the linear approximation with respect to the small dissipative param- 
eter 6 one can obtain an analytical expression for D(t). Keeping only first- 
order terms in the expressions for /x 2 , a 2 , a 3 and r„ arf and substituting them 
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into Eq. (157)) yields the diffusion coefficient: 

D = D °^( 1 + ^6»°&)> < 42 > 

where, as previously, 5' = 5 (2n(r)) 1 ^ 10 and c^o ~ 6.485. Note that due to the 
account of the third Sonine coefficient in /j, 2 and T Via d the linear approximation 
(I42p for D(t) differs from the previously obtained result 

The linear approximation ( 142)) is compared in Fig. [8] with the complete 
numerical solution. As it follows from the figure, after about ten collision 
per particle the two solution become practically indistinguishable. As it may 
be seen from the inset of Fig. [8], the impact of the third Sonine coefficient 
on the behavior of self-diffusion coefficient is rather small even for the large 
value of dissipation parameter 5 = 0.3. 

5. Conclusion 

We study evolution of a granular gas of viscoelastic particles in a homoge- 
neous cooling state. We use a new expression for the restitution coefficient e, 
which accounts for the delayed recovery of the particle material at a collision 
and allows to model collisions with much larger dissipation, as compared to 
previously available result for e. We analyze the velocity distribution func- 
tion and the self-diffusion coefficient. To describe the deviation of the velocity 
distribution function from the Maxwellian we use the Sonine polynomial ex- 
pansion. In contrast to the commonly used approximation, which neglects 
all terms in the Sonine expansion beyond the second one, we consider explic- 
itly the third Sonine coefficient. We detect a complicated evolution of this 
coefficient and observe that it is of the same order of magnitude, with respect 
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to the (small) dissipative parameter 5, as the second coefficient. This con- 



tradicts the existing hypothesis 



341 ]. that the subsequent Sonine coefficients 



d2, CL3 . . ., a,k are of an ascending order of some small parameter, character- 
izing particles inelasticity. Similarly as for the case of a constant restitution 
coefficient, we obtain an indication of divergence of the Sonine expansion for 
large dissipation, 5 > 0.3. For the asymptotic long-time behavior we derive 
analytical expressions for both Sonine coefficients, which agree well with the 
numerical data. 

Using the obtained third Sonine coefficient we compute the self-diffusion 
coefficient D and derive an analytical expression for D in the linear, with 
respect to 5, approximation. We show that the complete solution approaches 
the approximate analytical solution after a transient time of about ten colli- 
sions per particle. We observe, that in spite of the importance of CI3 for an 
accurate description of the velocity distribution function, its impact on D is 
rather small. 

We also study the evolution of the high-energy tail of the velocity distri- 
bution. Using the equation for the time-dependent slope of the tail and the 
obtained Sonine coefficients we find the amplitude of the tail and the thresh- 
old velocity, which demarcates the main part of the velocity distribution and 
the high-energy tail. We find the analytical expression for the asymptotic 
behavior of the threshold velocity, which agrees well with numerics, and ob- 
serve, that in a course of time it shifts to larger values; this implies that the 
high-velocity tail becomes less pronounced. 

Such behavior of the threshold velocity, of the Sonine coefficients and 
of the coefficient of self-diffusion, naturally, manifests that the properties of 
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the system tend to those of a gas of elastic particles; our theory quantifies 
evolution towards this limit. 

References 

[1] Brilliantov N V and Poschel T 2004 Kinetic theory of Granular Gases 
(Oxford University Press) 

[2] Goldhirsch I 2003 Annual Review of Fluid Mechanics 35 267 

[3] Hinrichsen H, Wolf D E 2004 The Physics of Granular Media Berlin: 
Wiley 

[4] Poschel T and Luding S 2001 Granular Gases 564 of Lecture Notes in 
Physics Berlin: Springer 

[5] Poschel T and Brilliantov N V 2003 Granular Gas Dynamics 624 of 
Lecture Notes in Physics Berlin: Springer 

[6] Goldsmit W 1960 The Theory and Physical Behavior of Colliding Solids 
(London: Arnold) 

[7] Bridges F G, Hatzes A and Lin D N C 1984 Nature 309 333 

[8] Kuwabara G, Kono K 1987 J. Appl. Phys. Part 1 26 1230 

[9] Tanaka T, Ishida T and Tsuji Y 1991 it Trans. Jap. Soc. Mech. Eng. 
57 456 

[10] Ramirez R, Poschel T, Brilliantov N V and Schwager T 1999 Phys. Rev. 
E 60 4465 

26 



[11] Brilliantov N V, Spahn F, Hertzsch J-M, Poeschel T 1996 Phys. Rev. E. 
53 5382 



[12 
[13 

[14 
[15 
[16 
[17 
[18 
[19 
[20 
[21 
[22 



Morgado W A M and Oppenheim I 1997 Phys. Rev. E 55 1940 

Brilliantov N V, Albers N, Spahn F and Poeschel T 2007 Phys. Rev. E 
76 051302 

Schwager T and Poschel T 1998 Phys. Rev. E 57 650 

Esipov S E and P6schel T 1997 J. Stat. Phys. 86 1385 

Goldshtein A and Shapiro M 1995 J. Fluid. Mech. 282 75 

van Noije T P C and Ernst M H 1998 Granular Matter 1 57 

Brilliantov N V and Poschel T 2006 Europhys. Lett. 74 424 

Brey J J, Cubero D and Ruiz-Montero M J 1996 Phys. Rev. E 54 3664 

Montanero J M and Santos A 1999 Granular Matter 2 53 

Nakanishi H 2003 Phys. Rev. E 67 010301R 

Noskowicz S H, Bar-Lev O, Serero D, Goldhirsch I 2007 Europhys. Lett. 
79 60001 



[23] Brilliantov N V and Pdschel T 2000 Phys. Rev. E 61 5573 
[24] Brilliantov N V and Pdschel T 2003 Phys. Rev. E 67 061304 
[25] Poschel T and Brilliantov N V 2003 in Q p.131 
[26] McNamara S 1993 Physics of Fluids A 5 3056 

27 



Goldhirsch I and Zanetti G 1993 Phys. Rev. Lett. 70 1619 

Brito R and Ernst M H 1998 Europhys. Lett. 43 497 

Brilliantov N V, Saluena C, Schwager T and Poeschel T 2004 Phys. Rev. 
Lett. 93 134301 

Poschel T, Brilliantov N V and Formella A 2006 Phys. Rev. E 74 041302 

Poschel T, Brilliantov N V and Formella A 2007 Int. J. Mod. Phys. C 
18 701 

Schwager T and P6schel T 2008 Phys. Rev. E 78 051304 

Bodrova A S and Brilliantov N V 2009 Moscow University Physics Bul- 
letin 64 131 

Huthmann M, Orza J A and Brito R 2000 Granular Matter 2 189 
Barrat A , Loreto V and Puglisi A 2004 Physica A 334 513 
Garzo V 2004 Physica A 343 105 



28 



